clear*
set maxvar 11000

global date _2025_10_28
cd "/Users/atri0032/Dropbox/World Bank/Temperature, poverty, and inequality/Estimations/"

******************************************************************************************
********************************** Figure 3: Heterogeneity analysis by income
******************************************************************************************

// Income classification
import excel using "${date}/_data/OGHIST.xlsx", clear firstrow sheet(list)

reshape long y, i(code) j(year) string
destring year, replace
rename y income_classification

save "${date}/_data/income_classification.dta", replace

// Low income
    // --- POOR 2.15 ---
    use "${date}/_data/spid_for_analysis_v2.dta", clear
	
	sort code year
	merge m:1 code year using "${date}/_data/income_classification.dta"
	drop if _merge==2
	drop _merge
	keep if inlist(income_classification,"L","LM")
		
	local temp_bin bin3_1 bin3_2 bin3_3 bin3_4 bin3_5 bin3_6 bin3_7  bin3_9 bin3_10 bin3_11 bin3_12 bin3_13
	local controls ln_GDPppp Unempilo agri manuf popgr 
	
    reghdfe poor215 `temp_bin' rainfall_sum humid `controls', absorb(geo_code2_new year) cluster(geo_code2_new) keepsingle	

	forval i=1(1)13 {
	cap gen estimate`i' = _b[bin3_`i']
	cap gen lb`i' = _b[bin3_`i'] - invttail(e(df_r),0.025)*_se[bin3_`i']
	cap gen ub`i' = _b[bin3_`i'] + invttail(e(df_r),0.025)*_se[bin3_`i']
	}

	gen estimate8=0
	gen lb8=0
	gen ub8=0

    gen n = _n
    keep if n == 1
    reshape long estimate lb ub, i(n) j(index)
    keep estimate lb ub index
    gen poor = 1

    tempfile poor215
    save `poor215'

    // --- Gini ---
    use "${date}/_data/spid_for_analysis_v2.dta", clear

	sort code year
	merge m:1 code year using "${date}/_data/income_classification.dta"
	drop if _merge==2
	drop _merge
	keep if inlist(income_classification,"L","LM")	
		
	local temp_bin bin3_1 bin3_2 bin3_3 bin3_4 bin3_5 bin3_6 bin3_7  bin3_9 bin3_10 bin3_11 bin3_12 bin3_13
	local controls ln_GDPppp Unempilo agri manuf popgr 

    reghdfe gini `temp_bin' rainfall_sum humid `controls', absorb(geo_code2_new year) cluster(geo_code2_new) keepsingle	

	forval i=1(1)13 {
	cap gen estimate`i' = _b[bin3_`i']
	cap gen lb`i' = _b[bin3_`i'] - invttail(e(df_r),0.025)*_se[bin3_`i']
	cap gen ub`i' = _b[bin3_`i'] + invttail(e(df_r),0.025)*_se[bin3_`i']
	}

	gen estimate8=0
	gen lb8=0
	gen ub8=0

    gen n = _n
    keep if n == 1
    reshape long estimate lb ub, i(n) j(index)
    keep estimate lb ub index
    gen poor = 2

    tempfile gini
    save `gini'	

    // --- Combine and plot ---
    use `poor215', clear
    append using `gini'

    replace index = index + 0.15 if poor == 2

    graph twoway ///
        (rcap lb ub index if poor == 1, pstyle(ci) lcolor(navy%80) msize(medium) lwidth(medium)) ///
        (scatter estimate index if poor == 1, mcolor(navy%80)) ///
        (rcap lb ub index if poor == 2, pstyle(ci) lcolor(maroon%80) msize(medium) lwidth(medium)) ///
        (scatter estimate index if poor == 2, mcolor(maroon%80)) ///	
        , scheme(white_tableau) ///
        xtitle("Temperature bin") ///
        xlabel(1 "<0" 2 "[0,3)" 3 "[3,6)" 4 "[6,9)" ///
	5 "[9,12)" 6 "[12,15)" 7 "[15,18)" 8 "[18,21)" 9 "[21,24)" 10 "[24,27)" 11 "[27,30)" ///
	12 "[30,33)" 13 "33+", nogrid) ///
        ylabel(, nogrid) ///
        ytitle("Effects of temperature") title("Panel A: Lower-income countries") ///
        yline(0, lcolor(black) lpattern(solid) lwidth(thin)) ///
        legend(row(1) order(1 "Headcount poverty" 3 "Inequality") pos(6)) name("gr1", replace)

    graph export "${date}/_figures/_fig_3_lowincome.pdf", as(pdf) replace	
	
	
	
// High income
    // --- POOR 2.15 ---
    use "${date}/_data/spid_for_analysis_v2.dta", clear
	
	sort code year
	merge m:1 code year using "${date}/_data/income_classification.dta"
	drop if _merge==2
	drop _merge
	keep if inlist(income_classification,"H","UM")
		
	local temp_bin bin3_1 bin3_2 bin3_3 bin3_4 bin3_5 bin3_6 bin3_7  bin3_9 bin3_10 bin3_11 bin3_12 bin3_13
	local controls ln_GDPppp Unempilo agri manuf popgr 

    reghdfe poor215 `temp_bin' rainfall_sum humid `controls', absorb(geo_code2_new year) cluster(geo_code2_new) keepsingle	

	forval i=1(1)13 {
	cap gen estimate`i' = _b[bin3_`i']
	cap gen lb`i' = _b[bin3_`i'] - invttail(e(df_r),0.025)*_se[bin3_`i']
	cap gen ub`i' = _b[bin3_`i'] + invttail(e(df_r),0.025)*_se[bin3_`i']
	}

	gen estimate8=0
	gen lb8=0
	gen ub8=0

    gen n = _n
    keep if n == 1
    reshape long estimate lb ub, i(n) j(index)
    keep estimate lb ub index
    gen poor = 1

    tempfile poor215
    save `poor215'

    // --- Gini ---
    use "${date}/_data/spid_for_analysis_v2.dta", clear

	sort code year
	merge m:1 code year using "${date}/_data/income_classification.dta"
	drop if _merge==2
	drop _merge
	keep if inlist(income_classification,"H","UM")	
		
	local temp_bin bin3_1 bin3_2 bin3_3 bin3_4 bin3_5 bin3_6 bin3_7  bin3_9 bin3_10 bin3_11 bin3_12 bin3_13
	local controls ln_GDPppp Unempilo agri manuf popgr 

    reghdfe gini `temp_bin' rainfall_sum humid `controls', absorb(geo_code2_new year) cluster(geo_code2_new) keepsingle	

	forval i=1(1)13 {
	cap gen estimate`i' = _b[bin3_`i']
	cap gen lb`i' = _b[bin3_`i'] - invttail(e(df_r),0.025)*_se[bin3_`i']
	cap gen ub`i' = _b[bin3_`i'] + invttail(e(df_r),0.025)*_se[bin3_`i']
	}

	gen estimate8=0
	gen lb8=0
	gen ub8=0

    gen n = _n
    keep if n == 1
    reshape long estimate lb ub, i(n) j(index)
    keep estimate lb ub index
    gen poor = 2

    tempfile gini
    save `gini'	
	
    // --- Combine and plot ---
    use `poor215', clear
    append using `gini'

    replace index = index + 0.15 if poor == 2

    graph twoway ///
        (rcap lb ub index if poor == 1, pstyle(ci) lcolor(navy%80) msize(medium) lwidth(medium)) ///
        (scatter estimate index if poor == 1, mcolor(navy%80)) ///
        (rcap lb ub index if poor == 2, pstyle(ci) lcolor(maroon%80) msize(medium) lwidth(medium)) ///
        (scatter estimate index if poor == 2, mcolor(maroon%80)) ///	
        , scheme(white_tableau) ///
        xtitle("Temperature bin") ///
        xlabel(1 "<0" 2 "[0,3)" 3 "[3,6)" 4 "[6,9)" ///
	5 "[9,12)" 6 "[12,15)" 7 "[15,18)" 8 "[18,21)" 9 "[21,24)" 10 "[24,27)" 11 "[27,30)" ///
	12 "[30,33)" 13 "33+", nogrid) ///
        ylabel(, nogrid) ///
        ytitle("Effects of temperature") title("Panel B: Higher-income countries") ///
        yline(0, lcolor(black) lpattern(solid) lwidth(thin)) ///
        legend(row(1) order(1 "Headcount poverty" 3 "Inequality") pos(6)) name("gr2", replace)

		
    graph export "${date}/_figures/_fig_3_highincome.pdf", as(pdf) replace	
	